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Abstract 

We investigate the gel formation from the equihbrium sol phase in a simple model 
that has the characteristics of (colloidal) gel-forming systems at a finite temperature. 
At low volume fraction and low temperatures, particles are linked by long-living 
bonds and form an open percolating network. By means of molecular dynamics 
simulations, we study the lifetime of bonds and nodes of the gel network in order 
to relate these quantities to the complex relaxation dynamics observed. 
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1 Introduction 



Gels are a special case of slow and disordered systems: Due to the low volume 
fraction, their slow dynamics are intimately related to the formation of an open 
structure (jH). The understanding of the connection between the structure of 
these systems and their complex dynamics is fundamental in the research on 
complex fluids, such as biological systems, polymeric materials, and colloids. 
The slow dynamics of gelling systems have some analogies with the dynamics 
of other disordered systems such as dense glasses, in that the majority of the 
particles becomes temporarily localized in space due to the structure forma- 
tion. On the other hand, due to the presence of the network, the gels show 
relaxation dynamics which strongly depend on the length scale/wave vector 
considered. At small wave vectors and close to the gelation threshold one 
typically observes very slow stretched-exponential and power law relaxations. 
Surprisingly, at higher wave vectors, recent experimental work on colloidal 
gels shows the onset of a compressed exponential relaxation (i.e. Kohlraush 
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exponent /3 > 1) This non-trivial behavior is the subject of an intense 
debate in the last years. In particular, it is still not clear how is it related to 
the structural features of the gel. 

Numerical simulation studies offer an effective tool in that they allow to fol- 
low the local relaxation processes and to understand their role in the complex 
dynamics. In spite of this, progress has also been hampered by the fact that 
there are so far very few microscopic realistic models for gels that effectively 
allow to study their dynamics by means of analytical methods or computer 
simulations. In colloidal systems the interparticle interactions can be suitably 
tuned: First, particles can be coated and stabilized leading to a hard sphere 
behavior. Then, an attractive depletion interaction can be induced by adding 
some non-adsorbing polymers. The range and strength of the potential are 
controlled by the size and concentration of the polymer respectively (jl). For 
volume fractions from 3% up to 40% and attractive interaction energies of the 
order of ksT confocal microscope images (0) show the formation of persis- 
tent structures, different at different volume fractions, and gel phases, whose 
mechanisms of formation are still poorly understood. It is important to note 
that here the underlying thermodynamics may significantly interplay and/or 
compete with gel formation via phase separation or microphase separation. 
As a consequence, coarsening or ordering processes will interfere with the gel 
dynamics. Whereas in the experiments the time scales typical of the micro 
or macro-phase separation are often much longer than the observation time 
scales, this is not the case in numerical studies using traditional models for 
colloidal suspensions, where the investigation of the gel dynamics has been 
severely hindered until very recently. 

In order to overcome these difficulties, we have recently proposed a new 
model (j^; 0) consisting of meso-particles which interact via a short range at- 
traction of the type already used to study colloidal systems. We have started 
from the idea that, in order to obtain a persistent open structure, a mech- 
anism competing with phase separation processes is needed. The novelty of 
our model is to include a directional effect in the interactions and therefore to 
facilitate the formation of the open structures found in gelling systems. These 
directional interactions are given by a combination of site-site interactions of 
the particles as well as a short ranged three-body interaction. Thus this new 
approach will promote the formation of large-scale disordered structures that 
are stable also at relatively high temperatures, without imposing a fixed local 
functionality of the meso-particle. It should also be noted that the presence 
of directional effective interactions, although not yet thoroughly investigated, 
is likely to play an important role in the physics of colloidal gels. This is sup- 
ported by the experimental observation of the local rigidity of the structures 
and of their low local connectivity: The open structures which spontaneously 
form at low volume fractions are characterized by aggregating units (particles 
or aggregates of particles) with a coordination number between 3 and 4 (0). 
In addition, surface inhomogeneities or sintering of the particles which are 
present in many cases can produce quite rigid links (jsj). 
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In this paper, we analyze the formation of the gel network in this model and 
focus on the time autocorrelation of bonds and network nodes. This study is 
essential to make the connection between the features of the network structure 
and the onset of the complex dynamics of the gel phase (@). 
In the following we briefly present the model and the numerical study in Sec. [2j 
We then describe the aggregation process and the structure of the gel phase in 
Sec. 131 In Sec. S] we study the time autocorrelation functions and the lifetime 
of bonds and of the nodes of the network structure. Finally, conclusions drawn 
from the results of the numerical study are contained in Sec. 



2 Model and numerical simulations 

The model investigated consists of identical particles of radius a, interacting 
via an effective potential Ves that is the sum of a two-and three body terms. 
The two body potential is itself the sum of a hard core like interaction that is 
given by a generalized Lennard- Jones potential, Vi^j{r), and a term V^p that 
depends on the relative orientation of the particles. For the radial term we 
have used VLj(r) = 23e[(o"/r)^^ — (a/r)^^] where the prefactors and exponents 
have been chosen in such a way to give a relatively narrow well of depth e. 
In the following we will measure length and energy such that a = 0.922 and 
e = 1, respectively, and time in units of ^Jma"^ /e, where m is the mass of a 
particle. 

As already discussed in the previous section, we have introduced a mechanism 
that favors the formation of an open network structure, in competition with 
the phase separation induced by purely radial interactions. This has allowed us 
to obtain an open persistent structure at temperatures and volume fractions 
where phase separation or microphase separation do not occur and therefore 
to study the dynamics of the gel formation without the complex interplay 
with the underlying thermondynamics. In practice, we have introduced a di- 
rectional interaction by decorating each particle with 12 points that form a 
rigid icosahedron inscribed in a sphere of diameter l.lcx. The potential 
between a particles #1 with a particle #2 is then set up in such a way that 
it is more favorable that the center of particle #1 approaches particle #2 in 
the direction of one of the points of the icosahedron that decorate particle 
^2. In addition we have also included an explicit (short range) three body 
potential V3 in the form of a gaussian in the angle 9 between three neighbor- 
ing particles and which makes that values of 9 smaller than 0.47r are unlikely. 
Although recently 3-body interactions have been directly measured amongst 
charged colloidal particles (jof), one could more generally see the particles in 
our model as the building blocks of the gel network (i.e. particles or aggre- 
gates of particles). In this terms, our model is aimed to describe the presence 
of directional effective interactions at some mesoscopic length scale, leading 
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to the formation of the open structures. More details on these potentials will 
be given elsewhere (jiol ). but here it should be noted that, whereas V3 is a real 
3-body term introducing an angular rigidity, is a geometric term which 
allows to add a soft-sphere repulsion depending on relative orientation. This 
term is not able alone to effectively limit the functionality of the particles at 
the volume fraction and temperatures considered here. The choice of consider- 
ing the effect of both these terms has been made in the spirit of investigating 
more deeply their relative contribution to the formation of the open network 
(and to its dynamics). For example, by taking only Vu -|- V3 an open struc- 
ture can certainly be obtained, but probably slightly different from the one 
discussed here. For the moment a complete understanding of the interplay of 
these two terms has not been reached yet. 

On the whole, in our model the (meso) particles can form directional bonds 
that favor the formation of an open network structure, without, however, im- 
posing a local symmetry or connectivity, in contrast to models that have been 



proposed before (I 111 : Il2l ). As a consequence, in our case there is an effective 



functionality of the particles, which is the result of the competition between 
entropic and internal energy contribution in the network formation. 
Using these interactions, we have done microcanonical simulations using con- 
strained molecular dynamics and a suitable combination of the algorithms 
RATTLE and SHAKE (1131 ) with a step size of 0.002. The number of particles 
was 8000 and the size of the simulation box L = 43.09 which gives a volume 
fraction of 0.05. (This corresponds to a particle density of 0.1) Before starting 
the production runs we carefully equilibrated the system by monitoring that 
the relevant time correlation functions have attained their asymptotic limit. 
The temperatures investigated were 5.0, 2.0, 1.0, 0.7, 0.5, 0.3, 0.2, 0.15, 0.1, 
0.09, 0.08, 0.07, 0.06, and 0.05. In order to improve the statistics of the results 
we have averaged them over five independent runs. 



3 Structure and cluster formation 



The change of the topology of the structure with decreasing T can be char- 
acterized by investigating the coordination number c(n), which is shown in 
Fig. [1] as a function of the inverse temperature. We define c{n) as the fraction 
of particles that have exactly n neighbors. (Two particles are considered to be 
neighbors if their distance is less that rmin=l.l, the location of the first mini- 
mum in the radial distribution function.) We see that at high temperatures a 
large majority of the particles are isolated {n = 0). By lowering the tempera- 
ture, the number of these free particles rapidly decreases and they practically 
disappear. The fraction of dimers or chain ends, n = 1, is around 30% at high 
T. With decreasing T such fraction increases, it attains a maximum at around 
T = 0.15, and then decreases quickly with decreasing T. At the same time the 
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Fig. 1. (Color online) T— dependence of c(n), the fraction of particles having coor- 
dination number n. 

number of particles that have exactly two nearest neighbors increases rapidly 
and these (local) configurations become by far the most prevalent ones at low 
T. Last not least also the number of particles with n = 3 neighbors increases 
quickly with decreasing T. From these curves we thus can conclude that with 
decreasing T the system forms an open network in which most particles form 
chains that meet at points with coordination number three and which are thus 
important for the mechanical properties of the structure. We have also shown 
elsewhere (jsl; 0) that whereas at high T the static structure factor S{q) of the 
system is relatively flat, i.e. the system has a structure similar to the one of 
a gas of free particles, around T = 0.1 S{q) starts to show a peak at go = 7.7, 
a wave-vector which corresponds to the distance between two nearest neigh- 
bor particles. In addition S{q) shows an increase at small q which indicates 
the formation of a disordered open network structure. This increase is rela- 
tively moderate and does not strongly depend on T, giving evidence that the 
system does not undergo a phase separation even at the lowest temperatures 
considered. Moreover, S{q) does not show any pronounced peak at these low 
wave-vectors and thus we can conclude that the network is disordered and 
does not have a well defined length scale. 

The appearance of the peak at go = 7.7 is due to the fact that at intermediate 
and low T the particles condense into clusters. By defining that a particle 
belongs to a cluster if its distance from at least one member of the cluster 
is less than rmin, we have studied n(s), the number of clusters that have ex- 
actly s particles (i5i). For T > 0.3 this distribution follows an exponential law. 
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corresponding to the random formation of transient clusters of non-bonded 
particles at low densities. Coherently with the data shown in FiglH at T = 0.1 
the shape of the distribution has strongly changed in that now the most prob- 
able clusters have size s = 2. At the same time n{s) has grown a tail at large 
s in that the largest clusters found have now O(IOO) particles. At lower tem- 
peratures, n{s) crosses over to a power law regime for high values of s, with 
a crossover point that moves to larger s with decreasing T. At T=0.06 this 
regime is apparently compatible with random percolation (jl4j ). At the lowest 
temperatures, T < 0.06, the distribution shows a gap at large s in that the 
system can form one big cluster that contains a substantial fraction of the 
particles in the system. Finally at T = 0.05 we have only very few particles 
that are members of small clusters, with an n{s) that is basically a constant, 
whereas the overwhelming majority (> 97%) of the particles belongs to one 
large percolating cluster. Hence, at very low temperatures it is energetically 
very unfavorable to have free particles and small free clusters. 
In Fig. [2] a typical configuration at T = 0.05 is shown: An open network in 
which most particles are forming chains (light spheres) that are connected in 
a relatively random way at points that have three or more neighbors (dark 
spheres) . 
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4 Dynamics: Lifetimes of bonds and nodes of the network 



Having characterized the structure of the system, we now turn our attention to 
its dynamical properties. We have studied the mean-squared displacement of 
a tagged particle, (r^(t)) = (|rj(t) — rj(0)R. A detailed discussion of the time 
dependence of can be found in Ref.(j6|). Using this function and the Ein- 

stein relation we obtained the diffusion constant D{T). Furthermore we have 
calculated the self-intermediate scattering function Fs{q,t) for wave- vector q, 
Fsiq,t) = N~'^ Z]jLi(exp[iq ■ (rj(t) — rj(0))]). The discussion of the quite com- 
plex q and T dependence of and Fs{q,t), characteristic for a gel- forming 
system, has also been presented elsewhere (l6|). From these functions we have 
calculated the characteristic relaxation times Ts{q,T) = J Fs{q,t)dt. 
Here we focus on the time autocorrelation functions of bonds and nodes form- 
ing the gel network. This analysis is essential to understand the role of different 
parts of the structure in the complex dynamics observed. As in Fig 12] we make 
the distinction of bonds that are connected to particles that have a coordina- 
tion number of three (referred to as "anchor particles") from the bonds that 
are connected to particles that have a coordination number of two ( "bridging 
particles"). For this we have determined Cb{t), the probability that a bond 
that exists at time zero is still present at time t, defined as 

^^^^^ ^ E., [(n.,(t)n,,(0)) - 



where nij{t) = 1 if particles i and j are linked at time t and nij{t) = 
otherwise. We have also calculated the corresponding time correlation function 
for the anchor points Csbit) defined as 

„ _ Ei l{n:dt)nu(0)) - {nx}''] 



where n^i{t) = 1 if particles i is an anchor point at time t, n^i{t) = otherwise. 
We have determined the corresponding average lifetime of bonds and anchor 
points as Tb(T) = J Ch{t)dt and r^bii) = J C3b(t)dt, respectively. In Fig. [3]these 
quantities are plotted in an Arrhenius plot together with the relaxation times 
Ts{q, T) obtained from the incoherent scattering functions Fs{q,t) for different 
wave- vectors q and the inverse of the diffusion constant D{T) ([5|). This figure 
well elucidates the role of bonds and nodes in the fundamental difference be- 
tween the high T and the low T regime of the relaxation dynamics: At low 
temperatures (i.e. T < 0.1) the bond lifetime sets the longest relaxation time 
scale in the system, and it becomes much longer than the simulation time 
window at the lowest T. That is, whereas at high temperatures the relaxation 
in the system is dominated by particle collision and diffusion, at the temper- 
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Fig. 3. (Color online) Arrhenius plot of the diffusion constant D and of the relaxation 
time Ts{q,T) as determined from the self-intermediate scattering function Fs{q,t). 
The solid line is a fit to the high T data for q = 3.0 of the form r = const. /v^- 
The open symbols correspond to different wave-vectors. The stars and filled circles 
are Tf, and r^b, respectively. 

atures where the aggregation starts to produce persistent structures the main 
mechanism for the relaxation is the bond breaking. 

From the figure, we recognize that at high T the relaxation time Ts{q, T) fol- 
lows a T~°'^— dependence for all g, which can easily be understood from the 
ballistic motion of the free particles and the small clusters. (Recall that in 
this T— range the distribution of the cluster size is independent of T ({sl).) For 
wave-vectors that are large, q > 2.0, the T— dependence of Ts{q,T) at low 
temperatures is relatively weak but still somewhat stronger than T~°'^. This 
can be understood by realizing that on the corresponding length scales the 
particles can still undergo an (almost) ballistic motion, despite the fact that 
they are, at low T, connected to other particles, since the whole local struc- 
ture is moving ballistically. This type of motion is no longer possible if one 
considers length scales that become comparable to the size of the mesh of the 
network which is around 10 and thus corresponds to a g smaller than 1.0. For 
q = 0.55 we find at low temperatures a T— dependence of that is close to 
an Arrhenius law whereas for smaller wave-vectors we find an even stronger 
T— dependence. On length scales that are thus comparable or larger than the 
typical mesh size of the network the T— dependence is thus very similar to the 
one characterizing the slow dynamics of dense glasses. 
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Also included in Fig. [3] is the inverse of the diffusion constant D, scaled by a 
factor A = 49.2 in order to make it coincide with the value of Ts{q = 0.14, T = 
1.0). We find that in the whole T— range investigated this quantity follows very 
closely the T— dependence of Ts{q, T) for small q, which is evidence that in this 
system the relaxation of the structure is closely linked to the diffusive motion 
of the particles, in contrast to the behavior found in dense glass-forming sys- 
tems in which dynamical heterogeneities make that the T— dependence of D 
is weaker than the one of Tg (llSl ). This behavior corresponds to the lowest 
temperature dynamical regime associated to the presence of the persistent 
network. 

Since here we would like to point out in particular the features of bond and 
nodes correlation at low T, we note that for T < 0.1, Tb shows an Arrhenius 
dependence on the temperature, in agreement with the picture mentioned 
above. As shown in the figure, the lifetime of the bonds is around a factor 
of 20 larger than the one of the anchor points, but the T — dependence is 
the same. This result is reasonable since the anchor particles experience on 
average more mechanical stress than the bridging particles and, due to their 
higher coordination, have less possibilities to yield to this stress. Hence it is 
more likely that their bonds are broken. Fig 13] also shows that the low T regime 
associated to the presence of the network, whose complex djTiamics have been 
studied in Ref.(0), corresponds to u being the longest time scale and r^b be- 
coming of the order of the longest relaxation times in the system. At these 
temperatures, there are practically no free particles (see Figs{T] and [2]) and 
more than 97% of them belong to the spanning network. Therefore one can 
envision the relaxation process in this system that the connecting branches 
of the network detach from the anchor particles, the branch reorients and at- 
taches itself to a different branch, thus creating a new anchor point. Of course, 
the possibility that a branch breaks at a bridging particle can not be neglected 
completely since there are significantly more bridging particles than anchor 
particles. Thus this type of motion will contribute to the relaxation dynamics 
as well. 

Let us now see how the time dependence of bond and nodes correlations at 
different temperature can be connected to the scenario just presented. In the 
main panel of Fig. HJ Cb(t) is plotted as a function of the time for different 
temperatures T. Although the long time decay of these functions is well de- 
scribed by a single exponential law at all the temperatures, we can actually 
distinguish two different regimes, one of high temperatures (T > 0.1) and the 
low temperature regime, corresponding to T < 0.1. At high temperature the 
particle collisions promote uncorrelated bond breaking or formation leading 
to a short time decay of bond correlation whose amount monotonically de- 
creases with temperature. At long times the bond breaking due to the energy 
activation process will eventually lead to the complete decay of correlations. 
At low temperatures instead, the energy activation process appears to be the 
only relevant process in the decay of bond correlation (due to the less and less 
important role of particle collisions). 
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Fig. 4. (Color online) Main panel: Bond correlation function Cfe(t) for differ- 
ent temperatures T. Inset: Time correlation function of the anchor points of the 
network, C3f,(t), for different temperatures T. 

In the inset of Fig. HI C^bit) is plotted as a function of the time for different 
temperatures. At the higher temperatures (T > 0.1) at this volume fraction 
particles of connectivity 3 are extremely rare, as already shown by Fig. [H 
As long as a spanning network is not formed, i.e. at temperatures T > 0.06, 
the long time decay of time correlation functions C^bit) follows a simple ex- 
ponential law, whose characteristic relaxation time increases with decreasing 
T. In addition to this regime, at temperatures when a persistent spanning 
network is present in the system, i.e. T < 0.055, the decay of C^bit) becomes 
also stretched, with a stretching exponent /3 decreasing with T (/3 ~ 0.55 at 
T = 0.05). This indicates that, once the network is formed and it is sufficiently 
persistent, the breaking and formation of the nodes will be associated not only 
to the overcoming of an activation energy but also to some heterogeneous and 
cooperative dynamic process (IgI: [lol). 

The scenario emerging from these data can be interestingly connected to the 
detailed analysis of the wave vector dependence of Ts{q, T) presented in Ref.(@). 
There we have discussed the striking differences between the high temperature 
and low temperature regime of the time correlations of particle displacement. 
At high temperatures, by decreasing the wave vector, Ta{q, T) smoothly crosses 
over from a ballistic to a diffusive type of behavior with an exponential de- 
cay of time correlations. At low temperatures, instead, once that a permanent 
spanning network is present, two clearly different regime, respectively at high 
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and low wave vectors, can be individuated. We have been able to show that the 
relaxation at high wave vectors is due to the fast cooperative motion of pieces 
of the gel structure, whereas at low wave vectors the overall rearrangements 
of the heterogeneous gel make the system relax via a stretched exponential 
decay of the time correlators. The coexistence of such diverse relaxation mech- 
anisms is apparently determined by the formation of the gel network (i.e. the 
onset of the elastic response of the system) and it is characterized by a typical 
crossover length which is of the order of the network mesh size. 



5 Conclusion 



We have analyzed the structural and dynamical features in the gel formation 
by molecular dynamics simulations of a new model for colloidal gels, based 
on directional effective interactions. We have focused on the change in the 
topology corresponding to the formation of the gel network and on the time 
correlation of bonds and nodes of the gel network. Our analysis shows that, 
once that bonds are practically permanent over the simulation time window, 
the onset of the persistent network corresponds to the onset of some hetero- 
geneous and complex dynamics, which is directly related to the anchor points 
of the network. It is evident that the anchor points are very important for 
the mechanical stability of the structure. On length scales of the order of 
the interparticle distance the structure is quite flexible whereas on the scale 
of the mesh size of the network, which is given by the typical distance be- 
tween the anchor points, the structure becomes somewhat more rigid due to 
the constrains imposed by the enhanced connectivity. However, even on that 
length scale the system is relatively soft as compared to the one found in dense 
glasses. With decreasing temperature the connecting chains will become stiffer 
and hence give rise to an increased effective interaction between the anchor 
points del;©. 
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